Debugging

Import Data

corrcounts_merge <- readRDS("~/VersionControl/senescence_benchmarking/Data/corrcounts_merge.rds")
metadata_merge <- readRDS("~/VersionControl/senescence_benchmarking/Data/metadata_merge.rds")
SenescenceSignatures <- readRDS("~/VersionControl/senescence_benchmarking/CommonFiles/SenescenceSignatures_divided_newCellAge.RDS")

library(markeR)
library(ggplot2)
library(ggpubr)
library(edgeR)

Loading required package: limma

?markeR

Scores

?CalculateScores

ℹ Rendering development documentation for "CalculateScores"

df_logmedian <- CalculateScores(data = corrcounts_merge, metadata = metadata_merge, method = "logmedian", gene_sets = SenescenceSignatures)

senescence_triggers_colors <- c(
  "none" = "#E57373",  # Soft red  
  "Radiation" = "#BDBDBD",  # Medium gray  
  "DNA damage" = "#64B5F6",  # Brighter blue  
  "Telomere shortening" = "#4FC3F7",  # Vivid sky blue  
  "DNA demethylation" = "#BA68C8",  # Rich lavender  
  "Oxidative stress" = "#FDD835",  # Strong yellow  
  "Conditioned Medium" = "#F2994A",  # Warm orange  
  "Oncogene" = "#81C784",  # Medium green  
  "Lipid Accumulation" = "#E57373",  # Coral  
  "Calcium influx" = "#26A69A",  # Deep teal  
  "Plasma membrane dysruption" = "#D32F2F",  # Strong salmon  
  "OSKM factors" = "#FFB74D",  # Bright peach  
  "YAP KO" = "#9575CD"  # Deep pastel purple  
)

cellTypes_colors <- c(
  "Fibroblast" = "#FF6961",   # Strong Pastel Red  
  "Keratinocyte" = "#FFB347", # Strong Pastel Orange  
  "Melanocyte" = "#FFD700",   # Strong Pastel Yellow  
  "Endothelial" = "#77DD77",  # Strong Pastel Green  
  "Neuronal" = "#779ECB",     # Strong Pastel Blue  
  "Mesenchymal" = "#C27BA0"   # Strong Pastel Purple  
)

cond_cohend <- list(A=c("Senescent"), # if no variable is defined, will be the first that appears in the ggplot
                    B=c("Proliferative","Quiescent"))

PlotScores(ResultsList = df_logmedian, ColorVariable = "CellType", GroupingVariable="Condition",  method ="logmedian", ColorValues = cellTypes_colors, ConnectGroups=TRUE, ncol = 6, nrow = 2, widthTitle=20, y_limits = NULL, legend_nrow = 2,xlab=NULL, cond_cohend = cond_cohend)

wrap_title_aux <- function(title, width = 30) {
  if (nchar(title) <= width) {
    return(title)  # No need to wrap if it fits
  }

  wrapped_title <- ""
  while (nchar(title) > width) {
    # Find positions of capital letters and symbols near the wrap point
    capital_pos <- gregexpr("[A-Z]", title)[[1]]
    symbol_pos <- gregexpr("(_|-|:)", title)[[1]]

    # Check for symbol breaks within the last few characters (width - 5 to width)
    valid_symbol_breaks <- symbol_pos[symbol_pos >= (width - 5) & symbol_pos <= width]

    if (length(valid_symbol_breaks) > 0) {
      # If a suitable symbol is found, break at the first valid symbol
      break_at <- valid_symbol_breaks[1]
    } else {
      # If no suitable symbol, look for capital letters within the same range
      valid_capital_breaks <- capital_pos[capital_pos >= (width - 5) & capital_pos <= width]

      if (length(valid_capital_breaks) > 0) {
        # If a capital letter is found, break just before the capital letter
        break_at <- valid_capital_breaks[1] - 1
      } else {
        # If no suitable symbol or capital letter, break at width
        break_at <- width
      }
    }

    # Append the wrapped line
    wrapped_title <- paste0(wrapped_title, substr(title, 1, break_at), "\n")

    # Update title with the remaining text after the break
    title <- substr(title, break_at + 1, nchar(title))
  }

  # Add the remaining part of the title
  wrapped_title <- paste0(wrapped_title, title)

  return(wrapped_title)
}

Debugging

Individual Genes

• PCA with genes from signatures only

Violin Expression Plots

IndividualGenes_Violins(data = corrcounts_merge, metadata = metadata_merge, genes = c("CDKN1A", "CDKN2A", "GLB1","TP53","CCL2"), GroupingVariable = "Condition", plot=T, ncol=NULL, nrow=2, divide="CellType", invert_divide=FALSE,ColorValues=senescence_triggers_colors, pointSize=2, ColorVariable="SenescentType", title="Senescence", widthTitle=16,y_limits = NULL,legend_nrow=4, xlab="Condition",colorlab="") 

Using gene as id variables

Correlation Heatmap

CorrelationHeatmap(data=corrcounts_merge, 
                   metadata = metadata_merge, 
                   genes=c("CDKN1A", "CDKN2A", "GLB1","TP53","CCL2"), 
                   separate.by = "Condition", 
                   method = "pearson",  
                   colorlist = list(low = "#3F4193", mid = "#F9F4AE", high = "#B44141"),
                   limits_colorscale = c(-1,0,1), 
                   widthTitle = 16, 
                   title = "test", 
                   cluster_rows = TRUE, 
                   cluster_columns = TRUE,  
                   detailedresults = FALSE, 
                   legend_position="right",
                   titlesize=20)

Warning: Heatmap/annotation names are duplicated: pearson's coefficient
Warning: Heatmap/annotation names are duplicated: pearson's coefficient, pearson's coefficient
Warning: `legend_height` you specified is too small, use the default minimal height.
Warning: `legend_height` you specified is too small, use the default minimal height.
Warning: `legend_height` you specified is too small, use the default minimal height.

Expression Heatmaps

ROC/AUC

Cohen’s d

PCA with genes from signature only

#' @importFrom edgeR DGEList
#' @importFrom stats prcomp
#' @import ggplot2
#' @importFrom ggpubr ggarrange

plotPCA <- function(data, metadata, genes=NULL, scale=FALSE, center=TRUE, PCs=list(c(1,2)), ColorVariable=NULL,ColorValues=NULL,pointSize=5,legend_nrow=2, legend_position=c("bottom","top","right","left"),ncol=NULL, nrow=NULL){
  
  legend_position <- match.arg(legend_position)
  
  if (is.null(genes)){
    genes <-  row.names(data)
  }
  
  data <- data[row.names(data) %in% genes, , drop=F]
  
  if (!nrow(data)>1) stop(paste0("Error: Number of genes should be >1; In your data you have only found the gene ",genes))
  
  # Ensure metadata matches sample order if provided
  if (!is.null(metadata)) {
    colnames(metadata)[1] <- "Sample"
    rownames(metadata) <- metadata$Sample
    metadata <- metadata[colnames(data), , drop = FALSE]
    y <- edgeR::DGEList( log2(data+1), samples= metadata)  
  } else {
    y <- edgeR::DGEList( log2(data+1))  
  }
  

  
  nPCs <- max(unlist(PCs)) # get the maximum number of PC based on the user's choice
    
  PCAdata <- stats::prcomp(t(y$counts), scale=scale, center=center)
  PCAcounts <- PCAdata$x
  PCAcounts <- as.data.frame(PCAcounts)
  
  if (nPCs > ncol(PCAcounts)) stop("Error: Number of genes too low for number of chosen PCs. Please reduce number of PCs.")
      
  PCAcounts <-  cbind(PCAcounts[,1:nPCs],y$samples) 
  
  pltList <- list()
  
  for (pc in PCs){
    pc <- unlist(pc)
    ev = PCAdata$sdev^2
    pc_x <- round(100*ev[pc[1]]/sum(ev),2) 
    pc_y <- round(100*ev[pc[2]]/sum(ev),2) 
    
    plt <- ggplot2::ggplot(PCAcounts, ggplot2::aes_string(y = paste0("PC",pc[1]), x =  paste0("PC",pc[2])))
    
    
    # Add jittered points, optionally colored by ColorVariable.Default: Brewer Pallette "Paired"
    if (!is.null(ColorVariable)) {
      plt <- plt + ggplot2::geom_point(ggplot2::aes_string(fill = ColorVariable), size = pointSize, alpha = 0.5, shape=21, color="black")
    } else {
      plt <- plt + ggplot2::geom_point(size = pointSize, alpha = 0.5, shape=21, color="black", fill="#D8D8D8")  
    }
    
    # If ColorValues is provided, use a manual color scale; otherwise, if ColorVariable is provided,
    # use a default brewer palette.
    if (!is.null(ColorValues)) {
      plt <- plt + ggplot2::scale_fill_manual(values = ColorValues)
    } else if (!is.null(ColorVariable)) {
      plt <- plt + ggplot2::scale_fill_brewer(palette = "Paired")
    } 
    
    # Add axis labels (including variance)
    
    xlab <- paste0("PC",pc[1],": ",pc_x,"% variance")
    ylab <- paste0("PC",pc[2],": ",pc_y,"% variance")
    titleplot <- paste0("PC",pc[1], " vs PC",pc[2])
    
    plt <- plt + ggplot2::labs(fill = "", x = xlab, y = ylab, title=titleplot)
    
    # Change theme
    plt <- plt +
      ggplot2::theme_bw()+
    ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust=1),
          plot.title = ggplot2::element_text(hjust = 0.5),
          legend.position = "bottom")
    
    # Adjust legend rows if legend_nrow is specified
    if (!is.null(legend_nrow)) {
      plt <- plt + ggplot2::guides(fill = ggplot2::guide_legend(nrow = legend_nrow, position = legend_position))
    }
    
    # Add reference lines
    plt <- plt +
      ggplot2::geom_vline(xintercept=0, linetype="dotted") + 
      ggplot2::geom_hline(yintercept=0, linetype="dotted")  
    
    
    pltList <- c(pltList, list(plt))
    
    
    
  }
  
  n <- length(pltList)
  
  if(n==1){
    plt <- pltList[[1]]
  } else {
    
    # Determine grid layout
    if (is.null(ncol) && is.null(nrow)) {
      
      ncol <- ceiling(sqrt(n))
      nrow <- ceiling(n / ncol)
      
    } else if (is.null(ncol)){
      
      ncol <- ceiling(n / nrow)
      
    } else if (is.null(nrow)){
      
      nrow <- ceiling(n / ncol)
      
    }
    
    if (!is.null(ColorVariable)) {
      plt <- ggpubr::ggarrange(plotlist = pltList, ncol = ncol, nrow = nrow, common.legend = TRUE, align = "h")
    } else {
      plt <- ggpubr::ggarrange(plotlist = pltList, ncol = ncol, nrow = nrow, align = "h")
    }
    
    
    
  }
  
  print(plt)
  
  invisible(list(plt = plt,
              data = PCAcounts))
  
}